home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / NS_ROOT / NS_ROOTS.PAS < prev    next >
Pascal/Delphi Source File  |  1992-05-11  |  12KB  |  429 lines
  1. { $A+,B-,D-,F+,I-,L+,N-,R-,S-,V-}
  2.  
  3. { $IFOPT N-}
  4.    { $E-}
  5. { $ELSE}
  6.    { $E+}
  7. { $ENDIF}
  8.  
  9.  
  10. Unit NS_Roots;
  11. {*****************************************************************************}
  12. {
  13. { TNumMethods
  14. {  +-TRoots
  15. {
  16. {*****************************************************************************}
  17.  
  18. Interface {===================================================================}
  19.    Uses
  20.        WinTypes,
  21.       WinProcs,
  22.       WinDos,
  23.       WObjects,
  24.       Strings,
  25.       StdHdr;
  26.  
  27.  
  28.     Type
  29.        TFloatingPoint = Real;          { Format for all floating point numbers }
  30.  
  31.        PNumMethods = ^TNumMethods; {+++++++++++++++++++++++++++++++++++++++++++++++}
  32.       TNumMethods = Object(TObject)
  33.           ErrorCode : Byte;                        { Holds most recient error code }
  34.  
  35.             Constructor Init;
  36.             Function    Min(X1, X2 : TFloatingPoint) : TFloatingPoint; Virtual;
  37.             Destructor     Done; Virtual;
  38.       End; {Object}
  39.  
  40.  
  41.         PRoots = ^TRoots; {+++++++++++++++++++++++++++++++++++++++++++++++++++++}
  42.       TRoots = Object(TNumMethods)
  43.           X1          : TFloatingPoint;          { Initial lower bracket value }
  44.          X2          : TFloatingPoint;          { Initial upper bracket value }
  45.          Tolerance   : TFloatingPoint;   { Refine root until + or - Tolerance }
  46.          Iter        : LongInt;            { Iterations performed to converge }
  47.          IterMax     : LongInt;           { Maximum number of iterations allowed (MUST be a positive number!) }
  48.          RootValue   : TFloatingPoint;         { Value at the calculated root }
  49.  
  50.           Constructor Init(vX1, vX2, vTolerance : TFloatingPoint; vIterMax : LongInt);
  51.             Function    BrentRoots : TFloatingPoint;     Virtual;
  52.             Function    BisectionRoots : TFloatingPoint; Virtual;
  53.             Function    NewtonRoots : TFloatingPoint;    Virtual;
  54.             Function    fx(X : TFloatingPoint) : TFloatingPoint;      Virtual;
  55.             Function    fxPrime(X : TFloatingPoint) : TFloatingPoint; Virtual;
  56.          Destructor  Done; Virtual;
  57.         Private
  58.           OrigX1      : TFloatingPoint;
  59.          OrigX2      : TFloatingPoint;
  60.  
  61.           Function    Bracketed(vX1, vX2: TFloatingPoint): BOOLEAN;
  62.       End; {Object}
  63.  
  64.  
  65.     Const
  66.          nearzero : TFloatingPoint = 1.0e-20;
  67.  
  68.  
  69. Implementation {==============================================================}
  70.  
  71.     {++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++}
  72.     {+++ TNumMethods ++++++++++++++++++++++++++++++++++++++++++ TNumMethods +++}
  73.     {++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++}
  74.  
  75.    Constructor TNumMethods.Init;
  76.        Begin
  77.           TObject.Init;
  78.       End;
  79.    {EndConstructor}
  80.  
  81.  
  82.  
  83.     Function TNumMethods.Min;
  84.    {**************************************************************************}
  85.    { Calculates the minimum of two numbers                                    }
  86.    {**************************************************************************}
  87.         BEGIN
  88.               IF (X1 < X2) THEN
  89.              Min := X1
  90.          ELSE
  91.              Min := X2;
  92.          {EndIf}
  93.         END;
  94.    {EndFunction}
  95.  
  96.     Destructor TNumMethods.Done;
  97.        Begin
  98.       End;
  99.    {End Destructor}
  100.  
  101.  
  102.  
  103.    {+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++}
  104.    {+++ TRoots +++++++++++++++++++++++++++++++++++++++++++++++++++ TRoots +++}
  105.    {+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++}
  106.  
  107.       Constructor TRoots.Init;
  108.    {**************************************************************************}
  109.        Begin
  110. {         NumMethods.Init;}
  111.             X1        := vX1;
  112.          OrigX1    := vX1;
  113.          X2        := vX2;
  114.          OrigX2    := vX2;
  115.          Tolerance := vTolerance;
  116.  
  117.          {--------------------------------------------------------------------}
  118.          { We don't want negative itterations!  If vIterMax is less than 1    }
  119.          { set IterMax to the largest value a LongInt will hold.              }
  120.          {--------------------------------------------------------------------}
  121.          If (vIterMax < 1) Then
  122.              IterMax := 2147483647
  123.          Else
  124.              IterMax := vIterMax;
  125.          {EndIf}
  126.  
  127.  
  128.       End;
  129.    {EndConstructor}
  130.  
  131.     Function TRoots.fx(X: TFloatingPoint): TFloatingPoint;
  132.    {**************************************************************************}
  133.         BEGIN
  134.           Abstract;
  135.         END;
  136.    {EndFunction}
  137.  
  138.     Function TRoots.fxPrime(X: TFloatingPoint): TFloatingPoint;
  139.    {**************************************************************************}
  140.         BEGIN
  141.           Abstract;
  142.         END;
  143.    {EndFunction}
  144.  
  145.     Function TRoots.BrentRoots;
  146.    {**************************************************************************}
  147.         CONST
  148.          Finished : Boolean = False;
  149.           FPP      = 1.0e-11; { Machine floating point precision }
  150.  
  151.         VAR
  152.           result, CC, DD, EE, fa, fb, fc, Tol1, PP, QQ, RR, SS, XM: TFloatingPoint;
  153.  
  154.         BEGIN
  155.          Iter      := 0;
  156.           X1        := OrigX1;
  157.          X2        := OrigX2;
  158.          ErrorCode := 0;
  159.          fa        := fx(X1);
  160.          fb        := fx(X2);
  161.  
  162.               IF (NOT(Bracketed(fa,fb))) THEN
  163.                 ErrorCode := 1
  164.               ELSE
  165.                   BEGIN
  166.                     fc := fb;
  167.  
  168.                     REPEAT
  169.                       IF (NOT(Bracketed(fc,fb))) THEN
  170.                       Begin
  171.                                  CC := X1;
  172.                         fc := fa;
  173.                         DD := X2 - X1;
  174.                         EE := DD;
  175.                      End;
  176.                   {EndIf}
  177.  
  178.                       IF (ABS(fc) < ABS(fb)) THEN
  179.                           BEGIN
  180.                         {---------------------------------------------------}
  181.                         { Rename X1, X2, CC and adjust bounding interval DD }
  182.                         {---------------------------------------------------}
  183.                                 X1 := X2;
  184.                         X2 := CC;
  185.                         CC := X1;
  186.                                 fa := fb;
  187.                         fb := fc;
  188.                         fc := fa;
  189.                           END;
  190.                   {EndIf}
  191.  
  192.                       Tol1 := 2.0 * FPP * ABS(X2) + 0.5 * Tolerance; { Convergence check }
  193.                       XM := 0.5 * (CC-X2);
  194.  
  195.                       IF ((ABS(XM) <= Tol1) OR (ABS(fa) < nearzero)) THEN
  196.                           BEGIN
  197.                                 result := X2;
  198.                                 Finished := TRUE;
  199.                                 RootValue := fx(result);
  200.                           END
  201.                       ELSE
  202.                           BEGIN
  203.                                 IF ((ABS(EE) >= Tol1) AND (ABS(fa) > ABS(fb))) THEN
  204.                                     BEGIN
  205.                                       SS := fb/ fa;
  206.  
  207.                                       IF (ABS(X1 - CC) < nearzero) THEN
  208.                                           BEGIN
  209.                                             PP := 2.0 * XM * SS;
  210.                                             QQ := 1.0 - SS;
  211.                                           END
  212.                                       ELSE
  213.                                           BEGIN
  214.                                             QQ := fa/fc;
  215.                                             RR := fb /fc;
  216.                                             PP := SS * (2.0 * XM * QQ * (QQ - RR) - (X2-X1) * (RR - 1.0));
  217.                                             QQ := (QQ - 1.0) * (RR - 1.0) * (SS - 1.0);
  218.                                           END;
  219.                               {EndIf}
  220.  
  221.                                       IF (PP > nearzero) THEN
  222.                                   QQ := -QQ;
  223.                               {EndIf}
  224.  
  225.                                         PP := ABS(PP);
  226.  
  227.                                       IF ((2.0 * PP) < Min(3.0*XM *QQ-ABS(Tol1 * QQ), ABS(EE * QQ))) THEN
  228.                                  Begin
  229.                                                EE := DD;
  230.                                     DD := PP/QQ;
  231.                                  End
  232.                                       ELSE
  233.                                   Begin
  234.                                                DD := XM;
  235.                                     EE := DD;
  236.                                  End
  237.                               {EndIf}
  238.                                     END
  239.                                 ELSE
  240.                                     BEGIN
  241.                                       DD := XM;
  242.                                       EE := DD;
  243.                                     END;
  244.                         {EndIf}
  245.  
  246.                                 X1 := X2;
  247.                                 fa := fb;
  248.  
  249.                                 IF (ABS(DD) > Tol1) THEN
  250.                                   X2 := X2 + DD
  251.                         ELSE
  252.                                   IF (xm > 0.0) THEN
  253.                                  X2 := X2 + ABS(Tol1)
  254.                                   ELSE
  255.                                  X2 := X2 - ABS(Tol1);
  256.                            {EndIf}
  257.                         {EndIf}
  258.  
  259.                                 fb := fx(X2);
  260.  
  261.                                 Inc(Iter);
  262.  
  263.                                 IF (Iter >= IterMax) THEN
  264.                               Begin
  265.                                      ErrorCode := 2;
  266.                                 Finished := True;
  267.                              End
  268.                            {EndIf}
  269.                           END;
  270.                   {EndIf}
  271.                     UNTIL (Finished);
  272.                   END;
  273.          {EndIf}
  274.  
  275.             BrentRoots := result;
  276.         END;
  277.    {EndMethod}
  278.  
  279.     Function TRoots.BisectionRoots;
  280.    {**************************************************************************}
  281.         CONST
  282.           FPP         = 1.0e-11;
  283.           nearzero = 1.0e-20;
  284.          Finished : Boolean = False;
  285.  
  286.         VAR
  287.           result     : TFloatingPoint;
  288.          fa, fb, fc : TFloatingPoint;
  289.          XC         : TFloatingPoint;
  290.  
  291.         BEGIN
  292.          Iter      := 0;
  293.           X1        := OrigX1;
  294.          X2        := OrigX2;
  295.               ErrorCode := 0;
  296.          fa        := fx(X1);
  297.          fb        := fx(X2);
  298.  
  299.               IF (NOT(Bracketed(fa, fb))) THEN
  300.                 ErrorCode := 1
  301.               ELSE
  302.                   BEGIN
  303.                     REPEAT
  304.                       XC := (X1 + X2) / 2.0;
  305.                   fc := fx(XC);
  306.  
  307.                       IF (Bracketed(fa, fc)) THEN
  308.                       X2 := XC
  309.                       ELSE
  310.                       Begin
  311.                                  X1 := XC;
  312.                         fa := fc;
  313.                      End;
  314.                   {EndIf}
  315.  
  316.                       IF (ABS(X2 - X1) <= Tolerance) THEN
  317.                           BEGIN
  318.                                 Finished  := TRUE;
  319.                                 result    := XC;
  320.                                 RootValue := fx(result);
  321.                           END
  322.                   Else
  323.                             Begin
  324.                               Inc(Iter);
  325.  
  326.                                 IF (Iter >= IterMax) THEN
  327.                               Begin
  328.                                    ErrorCode := 2;
  329.                                 Finished := True;
  330.                              End;
  331.                           {EndIf}
  332.                      End;
  333.                    {EndIf}
  334.                     UNTIL (Finished);
  335.                   END;
  336.          {EndIf}
  337.  
  338.               BisectionRoots := result;
  339.         END;
  340.    {EndMethod}
  341.  
  342.     Function TRoots.NewtonRoots;
  343.    {**************************************************************************}
  344.         CONST
  345.            Finished : BOOLEAN = False;
  346.           FPP        = 1.0e-11;
  347.             nearzero = 1.0e-20;
  348.  
  349.         VAR
  350.           result, lastX, fxCalc, FPX, XP : TFloatingPoint;
  351.  
  352.         BEGIN
  353.           X1        := OrigX1;
  354.          X2        := OrigX2;
  355.          ErrorCode := 0;
  356.          fxCalc    := fx(X1);
  357.          Iter      := 0;
  358.  
  359.               REPEAT
  360.                 IF (abs(X1) < nearzero) THEN
  361.                 Begin
  362.                     ErrorCode := 1;
  363.                   Finished := True;
  364.                End
  365.             ELSE
  366.                     BEGIN
  367.                       FPX := fxPrime(X1);
  368.  
  369.                       IF (abs(FPX) < nearzero) THEN
  370.                       Begin
  371.                           ErrorCode := 2;
  372.                         Finished := True;
  373.                      End
  374.                   ELSE
  375.                           BEGIN
  376.                                lastX := X1;
  377.                               X1 := X1 - fxCalc/FPX;
  378.                               fxCalc := fx(X1);
  379.  
  380.                               IF (ABS(X1 - lastX) <= Tolerance) THEN
  381.                                   BEGIN
  382.                                      Finished := TRUE;
  383.                                       result := X1;
  384.                                       RootValue := fx(X1);
  385.                                     END
  386.                         Else
  387.                                 Begin
  388.                                         Inc(Iter);
  389.  
  390.                                           IF (Iter >= IterMax) THEN
  391.                                         Begin
  392.                                              ErrorCode := 3;
  393.                                           Finished := True;
  394.                                        End
  395.                                      {EndIf}
  396.                                     End
  397.                         {EndIf}
  398.                            END
  399.                   {EndIf}
  400.                     END;
  401.             {EndIf}
  402.               UNTIL (Finished);
  403.  
  404.               NewtonRoots := result;
  405.         END;
  406.    {EndMethod}
  407.  
  408.     Function TRoots.Bracketed;
  409.    {**************************************************************************}
  410.    { Return True if vX1 * vX2 is negative.                                    }
  411.    {**************************************************************************}
  412.         BEGIN
  413.               IF ((vX1 > 0.0) AND (vX2 > 0.0)) OR ((vX1 < 0.0) AND (vX2 < 0.0)) THEN
  414.                 Bracketed := FALSE
  415.               ELSE
  416.                 Bracketed := TRUE;
  417.          {EndIf}
  418.         END;
  419.    {EndFunction}
  420.  
  421.    Destructor TRoots.Done;
  422.    {**************************************************************************}
  423.        Begin
  424.           TNumMethods.Done;
  425.       End;
  426.    {EndDestructor}
  427.  
  428.  
  429. End.